function [s_star, s_regions, s_hat, q_shat, q_sstar, q_F, dq_Fdd, ...
    T_star, T_smin, alpha_F_star, ...
    q_T, q_T_F, E_tax_F, E_tax_T_F, E_MC_tax_F, E_MC_tax_T_F, max_MC_tax,...
    Cdif, Cdif_i, Cdif_i_early, Cdif_i_late, Cdif_dep, Cdif_tax, ...
    MC_app0, MC_app0_i, MC_app0_i_early, MC_app0_i_late, MC_app0_dep, MC_app0_tax, ...
    MB_app0, MB_app0_i, MB_app0_i_early, MB_app0_i_late, MB_app0_dep, MB_app0_tax, ...
    dWdd_app0, dWdd_app0_i, dWdd_app0_i_early, dWdd_app0_i_late, dWdd_app0_dep, dWdd_app0_tax, ...
    W_app0, ind_max_W_app0] ...
    = fn_model_solution_app0(delta_grid, D0i_grid, R1, lambda, s_min, s_max, D0, prob, dist_s, dist_d, d_min, ...
    d_max, delta_grid_stp, chi, phi, gam, kappa, y1, y2, Y1tau)


s_star    = zeros(length(delta_grid),1);
s_regions = zeros(length(delta_grid),3);
q_sstar   = zeros(length(delta_grid),1);
q_F       = zeros(length(delta_grid),1);

T_star       = zeros(length(delta_grid),1);
T_smin       = zeros(length(delta_grid),1);
alpha_F_star = zeros(length(delta_grid),1);

q_T          = zeros(length(delta_grid),1);
q_T_F        = zeros(length(delta_grid),1);
E_tax_F      = zeros(length(delta_grid),1);
E_tax_T_F    = zeros(length(delta_grid),1);
E_MC_tax_F   = zeros(length(delta_grid),1);
E_MC_tax_T_F = zeros(length(delta_grid),1);
max_MC_tax   = zeros(length(delta_grid),1);

Cdif         = zeros(length(delta_grid),1);
Cdif_i       = zeros(length(delta_grid),length(D0i_grid));
Cdif_i_early = zeros(length(delta_grid),length(D0i_grid));
Cdif_i_late  = zeros(length(delta_grid),length(D0i_grid));
Cdif_dep     = zeros(length(delta_grid),1);
Cdif_tax     = zeros(length(delta_grid),1);


MC_app0         = zeros(length(delta_grid),1);
MC_app0_i       = zeros(length(delta_grid),length(D0i_grid));
MC_app0_i_early = zeros(length(delta_grid),length(D0i_grid));
MC_app0_i_late  = zeros(length(delta_grid),length(D0i_grid));
MC_app0_dep     = zeros(length(delta_grid),1);
MC_app0_tax     = zeros(length(delta_grid),1);

n_g = round(length(delta_grid)/4);

%%

s_hat  = fn_s_hat(R1,lambda,s_min,s_max,phi); % Computes s_hat, which is indepdent of delta
q_shat = dist_s.cdf(s_hat); % Fundamental failure
dq_Fdd = fn_dq(delta_grid,R1,D0,lambda,prob,dist_s,s_min,s_max,dist_d,d_min,d_max,phi,delta_grid_stp); % Computes partial of q wrt delta

parfor j = 1:length(delta_grid)
    delta_j = delta_grid(j);
    
    % Equilibrium
    s_star(j)           = fn_s_star(delta_j,R1,D0,lambda,s_min,s_max,dist_d,d_min,d_max,phi);
    q_sstar(j)          = dist_s.cdf(s_star(j));
    s_regions(j,:)      = [q_shat, q_sstar(j) - q_shat, 1 - q_sstar(j)]; %failure, multiple equil, non-failure
    q_F(j)              = q_shat + prob*(q_sstar(j) - q_shat); % failure probability
    
    T_star(j)           = fn_tax(delta_j,R1,D0,chi,s_star(j),phi,dist_d,d_min,d_max); % T(s^*)
    T_smin(j)           = fn_tax(delta_j,R1,D0,chi,s_min,phi,dist_d,d_min,d_max); % T(s_min)
    alpha_F_star(j)     = fn_alphaF(delta_j,R1,D0,chi,s_star(j),phi,dist_d,d_min,d_max);
    
    [q_T(j), q_T_F(j), E_tax_F(j), E_tax_T_F(j),...
        E_MC_tax_F(j),E_MC_tax_T_F(j),max_MC_tax(j)] = ...
        fn_E_tax(R1,delta_j,D0,chi,phi,s_min,s_max,dist_s,dist_d,d_min,d_max,prob,s_hat,s_star(j),kappa); % prob(T>0) and E[T]
    
    % Marginal welfare
    [Cdif(j),Cdif_i(j,:),Cdif_i_early(j,:),Cdif_i_late(j,:),Cdif_dep(j),Cdif_tax(j),...
        MC_app0(j),MC_app0_i(j,:),MC_app0_i_early(j,:),MC_app0_i_late(j,:),MC_app0_dep(j),MC_app0_tax(j)] = ...
        fn_dWdd_main_app0(R1,delta_j,D0,chi,prob,phi,dist_s,s_min,s_max,dist_d,d_min,d_max,kappa,...
        s_hat,s_star(j),lambda,y1,y2,Y1tau,gam,D0i_grid);
    
    if rem(j,n_g)==0; disp(j); end
    
end

%% Marginal Benefit

% MB Approx 0
MB_app0         = -dq_Fdd.*Cdif;
MB_app0_i       = -dq_Fdd*ones(1,length(D0i_grid)).*Cdif_i;
MB_app0_i_early = -dq_Fdd*ones(1,length(D0i_grid)).*Cdif_i_early;
MB_app0_i_late  = -dq_Fdd*ones(1,length(D0i_grid)).*Cdif_i_late;
MB_app0_dep     = -dq_Fdd.*Cdif_dep;
MB_app0_tax     = -dq_Fdd.*Cdif_tax;

%% dWdd

% omega = 1, m(j,s)=1
dWdd_app0         = MB_app0         + MC_app0;
dWdd_app0_i       = MB_app0_i       + MC_app0_i;
dWdd_app0_i_early = MB_app0_i_early + MC_app0_i_early;
dWdd_app0_i_late  = MB_app0_i_late  + MC_app0_i_late;
dWdd_app0_dep     = MB_app0_dep     + MC_app0_dep;
dWdd_app0_tax     = MB_app0_tax     + MC_app0_tax;

%% Optimum

W_app0 = integrate_W(dWdd_app0,delta_grid,delta_grid_stp);

[~,ind_max_W_app0] = max(W_app0);

end
